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Abstract 

We investigate the stability domains of ground states of generalized Hubbard mod- 
els with next-nearest neighbour interaction using the optimum groundstate approach. 
We focus on the ^-pairing state with momentum P — and the fully polarized ferro- 
magnetic state at half-filling. For these states exact lower bounds for the regions of 
stability are obtained in the form of inequalities between the interaction parameters. 
For the model with only nearest neighbour interaction we show that the bounds for 
the stability regions can be improved by considering larger clusters. Additional next- 
nearest neighbour interactions can lead to larger or smaller stability regions depending 
on the parameter values. 



*Work performed within the research program of the Sonderforschungsbereich 341 (Koln- 
Aachen-Julich) 
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1 Introduction 



Correlation effects are of great importance in condensed matter physics. Super- 
conductivity and ferromagnetism are two important phenomena which can arise 
in an interacting many-body system. Theoretical investigations usually begin 
with choosing a suitable Hamiltonian. In general this Hamiltonian is too com- 
plex and must be reduced to a reasonable model which gives only a simplified 
description of reality Such simplifications make it even more desirable to ob- 
tain exact results and compare these with experimental data. In addition, they 
can be used to check the results from computer simulations and approximative 
methods. 

The simplest model of correlated electrons was introduced independently by 
Hubbard, Gutzwiller and Kanamori in 1963 as an attempt to describe the 
effect of correlations for d-electrons in transition metals 0, 0, ||. This model 
consists of two terms, one describes discrete quantum mechanical motion of 
electrons (hopping) and the other one the on-site Coulomb interaction between 
electrons. Nevertheless, the Hubbard model is one of the most important models 
in theoretical physics and is believed to exhibit various phenomena including 
metal-insulator transition, ferromagnetism, antiferromagnetism and supercon- 
ductivity. In spite of its simplicity only a few exact results are known. For 
instance, Lieb and Wu solved the one-dimensional (D = 1) model by using 
Bethe-Ansatz-techmqae Q . The other class of exact solutions belongs to the 
limiting case D — oo, where a dynamical mean-field approximation becomes 
exact H |6|. However, the situation becomes much more complicated in the 
lower dimensional cases. 

In recent years a new, non-pertubative method was developed by Brandt and 
Giesekus 0. The main idea is to start with a well-known ground state and 
then construct a corresponding Hamiltonian in the form of a projection opera- 
tor. This approach permits to include a large class of interaction parameters. 
A generalization of this method was presented by Strack and Vollhardt 
^ U . Ovchinnikov improved some of the results obtained previously by using 
a different method Jl(J] (see also fll]]). His approach is based on Gerschgorin's 
theorem which gives a lower bound for the ground state energy of the Hamil- 
tonian and thus complements the usual variational principle which gives upper 
bounds. A much simpler and clearer method was used by de Boer and Schad- 
SCHNEIDER This method is called Optimum Groundstate Approach and 

was introduced by Klumper, Schadschneider and Zittartz for spin models 
[fl3[ . The basic idea is to diagonalize a specially chosen local Hamiltonian and to 
make all the local states which are needed for the construction of a given global 
ground state also local ground states by choosing the interaction parameters 
appropriately. This approach leads to some inequalities between the interac- 
tion parameters which represent the minimal stability region of the investigated 
ground state. Due to this restriction only a subspace of the parameter space 
can be examined. 

Using a larger local Hamiltonian enables in a natural way the inclusion of more 
interactions which determine the stability conditions. In general, one finds an 
extension of the stability domain of the ground state. Independently, Szabo 
took this into account and improved some results obtained previously |l4[ ] . Ad- 
ditionally, he examined the behaviour of the stability domain in the presence 
of next-nearest neighbour interaction parameters. For instance, in the case of 
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77-pairing state with momentum P = n he verified a shrinking of the stability 
region for a small ratio between nearest and next-nearest neighbour hopping. 
In contrast to his numerical approach we shall investigate various ground states 
using analytical calculations. 



A Hamiltonian of a many-body system on an arbitrary lattice but with homoge- 
neous a- nearest neighbour interaction can be split up into local Hamiltonians, 
e.g. H = ^2jj h\j. The minimal cluster □ consists of only two nearest (a = 1) 
neighbour lattice sites (ij) and the corresponding local Hamiltonian is called 
bond Hamiltonian. The largest cluster contains obviously all lattice sites and 
can be expressed by hjj = H. For small clusters □ the local Hamiltonian hu 
can be diagonalized exactly. This limits the tractable cluster size. By adding 
a trivial constant to the Hamiltonian H, which never changes the physics, one 
can achieve that the lowest eigenvalue eo of hu vanishes, i.e. eo = 0. In this 
case the lowest eigenvalue E of H is either positive (E > 0) or zero (Eo = 0), 
because hu is a positive-semidefinite operator and the sum of such operators 
is also positive semidefinite. In the special case Eq = a local condition for 
finding a ground state |\&o) exists: 



This equivalence can be understood by considering ^Zrj^ol^nl^o) = (^o\H\^o) 
= 0. Since the hu are positive-semidefinite, all (^ol^nl^o) must vanish, which 
in turn implies (|l]). In the case Eq = 0, the global ground state consists only 
of ground states of the local Hamiltonian and no excited local states are in- 
volved. A ground state of this type is called optimum ground state. To obtain 
such ground states for a given system one must perform two steps. First, the 
ground states of the local Hamiltonian must be determined. Then one has to 
check whether a global ground state can be formed using only these local ground 
states. 

3 The generalized Hubbard model 

The Hamiltonian of the generalized Hubbard-Model on a D-dimensional, hyper- 
cubic lattice with L sites and homogeneous a-nearest neighbour interaction can 
be split up into local Hamiltonians hffl . Due to homogeneity all local Hamil- 
tonians are equal and can be divided into two parts. The first part contains 
hopping and interaction terms: 
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(for all □). 



(1) 




tot / ^(Cjg-CjCT + Cj^C, 



) 




C£cr)(^"i, — a H~ f^j,— cr) 



<7 




(2) 



+ V a (ni - l){hj - 1), 
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where the pairs (ij) denote a-nearest neighbours, for instance nearest (a = 1) 
and next-nearest (a = 2) neighbours. The fermion operators c ia and Ci a create 
and annihilate electrons with spin a G {til} a t site i which is associated with 
the single tight-binding Wannier orbital. The corresponding number operators 
are hi„ = cj^c^ and hi = + hii- The SU(2) spin operators are given 
by Sf — (h^ + h i i)/2, — ct.Cjf and — cLc,-^. The physical nature 
of the various terms is as follows: The first term (t) is the usual hopping of 
fcrmions on a lattice. The next two terms, bond-charge interaction (A) and pair- 
hopping (Y ) , were studied in relation with superconductivity [[l5[ |f6[ |l7| . The 
fourth term is an anisotropic Heisenberg term with a AAZ-type spin interaction 
given by the exchange constants J xy and J z . The last term (V) is known as 
the a-nearest neighbour Coulomb interaction. Estimates for the values of the 
couplings (for metals) for example can already be found in Hubbard's original 
paper [0. 

The second term contains only on-site interactions 0,-j = Oi + Oj with 

O t = | (ni T - l/2)(n iX - 1/2) + |rW. (3) 

Here U is the on-site Coulomb interaction, fi the chemical potential and Z the 
coordination number of nearest neighbour sites on the £>-dimensional hypercu- 
bic lattice. 

A local Hamiltonian hjj can be divided into bond Hamiltonians such that a com- 



parison with the results obtained in |12 is possible. We restrict our extension 



to cluster sizes A(D) = {3,4} and call the corresponding local Hamiltonians 
3- and ^-site Hamiltonian. In this case only nearest (a = 1) and next-nearest 
(a = 2) neighbour interactions exist on the square lattice and therefore: 

Aa:=| E {^ + O i0 ) + £ (4) 

The factor F := 2(D — 1) for D > 1 is only necessary in order to compare results 
of different clusters without rescaling the coupling constants (due to multiple 
counts of bonds). Fig. |l| shows the covering of 5-site clusters on a square lattice. 
The entries of a local state are described by £ e {0, |, I, 2}, where '0' denotes 
an empty site, a £ {T, 1} a site occupied by one electron with spin a and '2' a 
doubly occupied site. The local state of cluster size N(\3) is a tensor product 
166 • ■ ■ 6v(D)) = 16) <S> 16) ® • • • ® I6v(d)). Together one gets local states 

and a x 4 Ar ( n )-matrix which represents the local Hamiltonian. Although 

this matrix might be very large, the number of zero elements is still a large 
number. The use of symmetries makes the problem more tractable. One of the 
simplest symmetries is associated with the conservation of the total number of 
electrons. One has to consider only subspaces corresponding to a fixed number 
of electrons, i.e. one gets a block diagonal matrix. Another useful condition 
which we shall frequently impose is X = t. It leads to the preservation of the 
number of doubly occupied sites (see e.g. H [18)) and thus some of the block 
matrices split into smaller ones. Table @ summarizes the results for the 5-site 
Hamiltonian with corresponding 64 x 64-matrix: 
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X ^ t 




X = t 




number 


block size 


number 


block size 


4 


1 x 1 


4 


1 x 1 


8 


3x3 


12 


3x3 


4 


9x9 


4 


6x6 



Table 1: Block sizes and number of blocks for the 3-site Hamiltonian. 



However, the determination of algebraic eigenvalues of a characteristic polyno- 
mial p(X) = det(M — AI) is limited, i.e. only polynomials up to fourth degree 
can be solved in closed form. In the case X = t it is possible to find a conve- 
nient transformation with a corresponding matrix T. After the transformation 
T _1 MT = M the four 6 x 6-matrices decay into blocks of size 3 and all eigenval- 
ues can be obtained in closed form. This is the main reason why we concentrate 
here on the case X = t. 



4 Results 

We shall restrict our discussion in the following to two physically interesting 
classes of states. The first class are the 77-pairing states with momentum P which 
show off-diagonal long-range order (ODLRO) and are thus superconducting pg[ . 
The second one is the fully polarized ferromagnetic state at half-filling. We 
determine under which circumstances these states are optimum ground states 
of the generalized Hubbard model. We shall mainly consider the square lattice 
(D = 2)._ 

The definition of an 77-pairing state with momentum P is given by the expression 

\v) = (v i p) 10) with vp = J2 eiPid ]Av ( 5 ) 

where Af is an integer which is related to the particle number N through Af = 
N/2. Since we would like the 77-pairing state to be the ground state of the global 
Hamiltonian it is informative to determine the commutator [H, 77J,]. A long, but 
straightforward calculation yields: 

2 

[H, tjU = 2(X Q - to) £ (e P] + e lFfe ) (cj^ + 

<*=1 (jk) a 



<J*>c 



(let iPj , r iPk\ , -t 

+ ( "y e ~ v <* e ) m - l ) c kAi 



(Jk) a v / 

(6) 
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Using (Q) one finds the conditions under which the r\— pairing states (||) are 
eigenstates of H. For the momenta P € {0, ir} we have the following constraints 
on the interaction constants: 



P = 


P = ir 


X a — t a 

Y a = 2V a 


Y a = (-1)<W Q 



Note that for momentum P = n no conditions relating ti to Xi exist. In the 
following we shall only consider the P — case. Other values of the momenta 
can be treated similarly. An investigation of the properties of these states can 
be found e.g. in fl8fl . 

In order to make the 77-pairing state (with momentum P = 0) an optimum 
ground state we first observe that |r;) can be built completely from the local 
3-site states 000), |222), 002) + |200) + |020) and |022) + |220) + |202) and 
analogous ^-site states. Without next-nearest neighbour interactions all local 
states have the same local energy eo = U/(2Z) + V if we set /i = 0. Demanding 
that eo is the local ground state energy and hence all other local energies must 
be larger, one obtains the following inequalities: 



^<0, J < min {»<">, ®< n \ »< n \ *< n >,...} (n=3,4) 
fg(s,4) ._ _ 2 | t | _ 2V 
4 

(g(S,4) . = _y + 1 ^_JZ „ ^/(J*)2+8(J»V)2) 



These inequalities represent the stability regions for the 77-pairing state with 
momentum P = 0. The selected^] bounds 58^ belong to the 5-site and/or 4~ 
site case and can be distinguished by the upper index n. There are seven bounds 
for n — 3 and more than 50 for n = 4 ■ The first two bounds were also obtained 
by using the bond-diagonalization . The other bounds are new and indicate 
an improvement of the stability region. 

It is possible to investigate all cuts of the parameter space (t,U,V, J z , J xy ), 
but we shall concentrate only on some of them. For all cuts we took realistic 
parameter values satisfying U>V>t>J z > J xy . 

In the J z — J xy cut of the parameter space (Fig. ||) the inner triangle corresponds 
to the stability region of the ^-pairing state with momentum P = obtained by 
bond-diagonalization. The enlargement corresponds to results achieved by the 
5-site Hamiltonian (including purely nearest neighbour interaction). The ^-site 
Hamiltonian yields no further improvement of the bounds. However, it is not 
clear that the ^-pairing state is not a ground state outside of those bounds since 
larger cluster sizes might yield a further enlargement of the stability region. In 

^^We have listed in (B) only those bounds which are relevant for figures ^ and j^. A complete 
list of the bounds 58^. S "can be found in the appendix. The bounds can be found in |2o| . 
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contrast to the last figure the following U — V cut (Fig. ||) displays also an 
enhancement achieved by ^-site diagonalization. 

The inclusion of next-nearest neighbour interactions modifies the local ground 
state energy eo = U / (2Z) + V1 + V2 and hence the constraints concerning the 
interaction parameters: 



Vi<-4V 2 , - <min{93i, Q5 2 , 58s, !B 4 ,...} 

I81 := -Vi - v 2 + ~ (J? + JJ) 
<B 2 := 2 f-ta -V1-V2- \J ' (t 2 + V 2 f + t{ 

33 3 := g (-8 (Vi + Va) - Ji + 2Jt v - \J ( Jf - 2 Jf + 2 J^) 2 + 8( J* B ) 2 

234 := i 2|«a| - 5 (Vi + V a ) + \ (J! - 3J|) + \ (Jf s + 3J|")) 

- -W( -8 I* 2 I + 4 (^1 + V2) - + 3J| + - 6JD 2 + 192t? 



These bounds belong to the 5-site diagonalization results since larger clusters 
cannot be diagonalized analytically in closed form. In order to be close to real 
systems we take smaller next-nearest neighbour parameters than corresponding 
nearest neighbour ones and express this through ratios rp :— P1/P2 with P a £ 
{t a , V a , J*, Ja V }- The ratios depend on the material and hence can be very 
different. Since the 77-pairing states with momentum P = consist of electron 
pairs it is interesting to consider the Y\ — Y 2 cut (Fig. |]). This cut represents 
the behaviour of the stability region for different on-site Coulomb interaction 
parameters U. All other parameter pairs have the same ratio rp = 3. Note 
that on the square lattice the numbers of nearest and next-nearest neighbours 
are exactly the same. One important observation is that the two parameters Yi 
and Y2 stabilize the ground state with increasing Coulomb repulsion (U > 0) 
which try to seperate the electron pairs. 

The fully polarized ferromagnetic state is a simple tensor product of local states 

L 

|F) = IK T I°>> w 

1=1 

where each lattice site is occupied by an electron with spin a =\ (at half- 
fillingQ). In contrast to the 77-pairing state this state is already an eigenstate of 
our Hamiltonian and hence there are no restrictions concerning the parameters. 
Nevertheless, we concentrate on the special case^] X = t. If we want \F) to 
become an optimum ground state we have to make \aoa) or \aoao) the local 
ground state. Without next-nearest neighbour interaction the corresponding 
local energy eo = —U/(2Z) + J z /4 — 2[i/Z is a lower bound for the other local 

2 Away from half-filling the ferromagnetic state is no optimum ground state. 
3 For X ^ t one has to rely on numerical methods. Preliminary results show a behaviour 
similar to the case X = t. 



(8) 
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energies leading to the inequalities: 



u 



> max 



03<"\ OS 



(n) 



} (n = 3,4) 



SB* 



(3,4) 



= 2|*| + T + 



03 



{a,4) 



z 

2\n\ 



(3) 



03 



£1 

~~2 



m 
z 



v • 



+ 2Y 2 



V + — + VV 2 + 2F 2 



(10) 



The first two boundaries are the same as those of the bond-results |12] ], the 
last ones are new. Hence an improvement of the ground state region might be 
obtained by considering the last two bounds. But the number of all possible 
two dimensional cuts is still very large. 

The bounds can be further improved using the following argument fl2]| . With a 
fixed particle number N a state is a ground state of H but also a ground state 
of H + /iN. In this situation one can regard the bounds as a function of \x and 
try to find the value of fj, which optimizes these bounds. For instance, if we have 
inequalities like a > b + fi and a > c — /z then the best value is /i = (c — b) /2 and 

(3) (4) 

thus a > (b + c)/2. In our case we get /i = and therefore 5$3 = QSjj . This 
leads to the cognition that only the 5-site diagonalization is necessary, and the 
inclusion of four local lattice sites does not improve the stability region. This 
result is shown in po| for various two dimensional cuts of the parameter space. 
With next-nearest neighbour interaction parameters we get the modified local 
energy eo = —U/(2Z) + Jf/4 + Jf/4 with the corresponding new inequalities: 



Ji < 
u 



-2(J| + Jf) 



> max {93 1, 03 2 , 93= 



(JiY 

93 4 ,.. 



(Jf Y < -2 (J! 



} 



OSi := -Vi - V 2 + - (Jf + JS) 

03 2 := 2t 2 + i (Ji + Ji) + i^/(4t 2 + J|) 2 + 16t? 



03:; 



-Vi - -Y a 



Ji 



+ ^\J(y 1 + Y 2 y + 2Y 2 - Wi (Vi + Y2-V2) 



03 4 := 2f 2 + Jr + J| - ± ( + J 2 I! 



+ ^(«2 + JTY + (J( - JD (JI - Ji y + 2J2 v + st 2 ) + mtj 



(11) 



The pure Hubbard model H — H(t±, U) exhibits no ferromagnetic ground state 
at half-filling. Only for some special cases like the Nagaoka case the ex- 
istence of a fully polarized ferromagnetic ground state can be proven. An ex- 
tension to the case of the generalized Hubbard model can be found in p2|]. 



The influence of long range hopping t a on ferromagnetism (at half-filling) was 
investigated e.g. by FARKASOVSKY ^3|. The results show a suppression of fer- 
romagnetism with increasing a. In the U — t\ cut (Fig. ^ we considered the 
behaviour of the stability region for different ti values. But in contrast to |^3| 
all other type of couplings were not turned off. The inclusion of the next-nearest 
neighbour hopping shows a reduction of the stability region for the ferromag- 
netic ground state in agreement with the results of [ [23| . 



5 Conclusions 

We have presented exact results for stability regions of two physically inter- 
esting ground states of the generalized Hubbard model with nearest and next- 
nearest neighbour interaction using the optimum ground state approach. First 
we looked at the ^-pairing state with momentum P = and then at the fully 
polarized ferromagnetic state at half-filling. We have studied the behaviour of 
the stability domains of these two states with increasing cluster size. But due 
to difficulties that emerge for analytical diagonalization, we have limited our 
analytical calculations to two cluster sizes, i.e. N(D) = {3,4}. These cluster 
Hamiltonians were divided into bond Hamiltonians so that a comparison with 
results obtained by bond diagonalization Q was possible. The new boundaries 
which exhibit an improvement were illustrated graphically in some chosen cuts 
of the parameter space. Without next-nearest neighbour interactions all cuts 
show an enlargement of the stability domains obtained by 5-site diagonalization. 
The extension to four lattice sites only have led to an amplification in the U — V 
cut, and all other cuts indicated a fast convergence of the stability regions. We 
expect that any further improvement is rapidly decreasing with increasing clus- 
ter size. Another aim of this work was to determine the effects of next-nearest 
neighbour interactions on the stability domains. However, we restricted our 
investigation to the 5-site case only. The illustration of these bounds indicates 
that the stability conditions are strongly dependent on the next-nearest neigh- 
bour parameters. For instance, in the ferromagnetic case one finds a reduction 
of the stability domain (U — t\ cut) in the presence of next-nearest neighbour 
hopping t2- In the case of the ^-pairing state with momentum P = we consid- 
ered the Y\ — Y% cut for different values of the on-site Coulomb interaction U. 
With increasing Coulomb repulsion (U > 0) which tries to separate the electron 
pairs we observed a stabilization of the domain because of large negative Y\ 
values. 

In summary, we have shown that a cluster of three and four lattice sites can 
be treated (with some restrictions) analytically whereby every two-dimensional 
cut of the whole parameter space can immediately be examined and does not 
require long numerical calculations. Although we restricted ourselves mainly to 
the case X = t, where the local Hamiltonian decays into small blocks which can 
be diagonalized in closed form, we like to stress that the optimum ground state 
approach can be used to treat the general case X ^ t as well. Here one can 
either diagonalize the larger block matrices numerically |l4| or combine the op- 
timum ground state method with the Gerschgorin approach of pOj . Instead of 
determining the eigenvalues of the larger matrices exactly one can obtain lower 
bounds in closed form by using Gerschgorin's theorem. Although these bounds 
will in general not be the best possible ones they still yield exact stability regions 
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for the state under consideration. The results give valuable information about 
the phase diagram. In Fig. |^ the stability regions of the two states investigated 
here are shown for a generalized Hubbard model with only nearest neighbour 
interactions. For the parameter values chosen one already knows a considerable 
part of the phase diagram. These results might serve as a guidance for further 
computer simulations or exact diagonalization studies. Apart from this aspect, 
the extension enables some interesting new investigations due to the inclusion 
of more correlations. 
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7 Appendix 

To obtain an optimum ground state for a given system one has to determine 
all local eigenvalues. For instance, in the 3-site case (with /i = 0, X a = t a and 
next-nearest neighbour interaction) the eigenvalues have the form: 





2Z 4 2 



e7,8 



U 



_i 

8 



-xy 



1 



yVf - 2 Jf + 2^ + 8(^)2, 



2Z 



+ 



4 



8 




ei3,i4 



U_ ta 
%Z T 2 



•II 
16 




8 8 4 
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ei5,i6 := 20 cos 
ei7,i8 : = 20 cos 
eig,20 := 20 cos 



^±2t 2T Y 2 -^-V 2 



9 + 2tt 
3 

6 + a-k 



i 

3 

-(- 

3 \4Z 



^±2t 2T Y 2 -^-V 2 



± 2t 2 =f y 2 



Vi 



^2 



The expressions of the functions O and 0, which include the interaction terms, 
are to large and hence are omitted. 

The bounds are derived from the inequalities e < (for all i). For the ry-paring 



state with P = one has e 



(r>) 



e\. It becomes an optimum ground state if the 



conditions V < and U/Z < min{*Bi, . . . , Q5 7 } are satisfied, where 



<8 



(3 



25 



<8 



(3 



(3 



23 



23< 5 



2 1*| - 2V, 



= -V 

1 

3 
1 

6 
1 

6 



~8 



■y/(J*) 2 + 8(J*v) 2 , 



J z 



\J xy \ 



- - J z - 2\j*y\y + 192* 2 ) 



2 4 

(-8V + J z ), 
(-8V-J' + 2\J xy \), 
I (-5V - J -- \^(AV + J z f + I92t^ 



Due to the complexity of the last six eigenvalues (ei5 - e 2 o) the corresponding 
bounds do not exist in "closed form", but numerical investigations show that 
they are irrelevant. Including also next-nearest neighbour terms one gets: 



Si 

<B 7 



= -v 1 -v 2 + -u z l +j 2 z ) 1 



2 (-* 2 -v 1 -v 2 - ^(t 2 + v 2 y + t\ 

-8(Vi + K 2 ) - Jf + - 



2 J| + 2J 2 ry ) 2 + 8(Jf ) 2 



2*2 

3 



1 



|(^i + ^ - ^(jr + 3JD + -(j^ + 3j 2 ^) 



3 



12 

^ - + V 2 ) - l(Jf + 3J|) - i(J^ + 3JD 



1 



(8t 2 + 4(Vi + V2) - Jf + 3J| - + GJ* 2 ') 2 + 192* 2 



12 

4* 2 - 2Vi - 4V 2 , 

-Vi-v 2 - J(J 2 z + 2J 2 cy ), 
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®s = -k-8t 2 + 8Vi + 8V 2 - Jf), 
o 

®9 = -k8|t 2 | + 8Vi + 8V 2 + Jf +2|Jf|), 
o 



- ^ \/ (_8 * 2 + 404 + V^) + Jf - 3J|) 2 + 192tf. 

The ferromagnetic state has = e 2 and the corresponding bounds can be 
derived in an analogous way. 

The number of eigenvalues and bounds in the case n = 4 is too large to be listed 
here. They can be found in EG]. 
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Figure 1 : Covering of 3-site clusters on a square lattice. In order to obtain the full 
lattice including next-nearest neighbour interaction (diagonal bonds) one has to cover 
the lattice with two different 3-site clusters, namely triangle of type 'lij' (upper part) 
and 'ijk' (lower part). 
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Figure 2: Stability region of the n-pairing state with momentum P = in the J z — J xy 
cut in units of \t\ for different cluster sizes. 
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Figure 3: A U — V cut for the rj-pairing state with momentum P = in units 



15 




Figure 4: A Yi — Y2 cui /or tfte bounds of the r/-pairing state with momentum P — 
twft non-zero next-nearest neighbour interactions. The stability domains are shown 
for different values of U . 
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Figure 5: Bounds for the stability region of the ferromagnetic state in the U — ti plane 
for different values of the next-nearest neighbour hopping ti . 
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